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Abstract 

We study the spinodal instabilities of asymmetric nuclear matter at finite temperature within 
the microscopic Brueckner-Hartree-Fock (BHF) approximation using the realistic Argonne V18 
nucleon-nucleon potential plus a three-body force of Urbana type. Our results are compared with 
those obtained with the Skyrme force SLy230a and the relativistic mean field models NL3 and 
TW. We find that BHF predicts a larger spinodal region. This result is a direct consequence of the 
fact that our Brueckner calculation predicts a larger critical temperature and saturation density 
of symmetric nuclear matter than the Skyrme and relativistic mean field ones. We find that the 
instability is always dominated by total density fluctuations, in agreement with previous results 
of other authors. We study also the restoration of the isospin symmetry in the liquid phase, i.e., 
the so-called isospin distillation or fragmentation effect, finding that its efficiency increases with 
increasing proton fraction and decreases as temperature and density increase. In general, we find 
that the Brueckner results are comparable to those obtained with the Skyrme and the relativistic 
mean field models, although the restoration of isospin symmetry is not so efficient in this case. 

PACS numbers: 21.65.-f, 21.30.-x, 25.70.Pq, 24.60.Ky 
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The nature of the nucleon-nucleon interaction gives rise to a nuclear matter equation of 
state of Van der Waals type and, therefore, one expects a liquid-gas like phase transition 
to take place l|. This transition can be explored in heavy-ion collisions; in particular the 
experimental results related to the multi-fragmentation phenomenon, where highly excited 
composed nuclei are formed in a gas of evaporated particles, can be interpreted as the co- 
existence of a liquid and a gas phase {2], 3], (4, 0, [(J. Since nucleons can be either neutrons or 
protons, one should consider nuclear matter as composed of two different fluids. Therefore 
the phase transition may lead to richer phase diagrams in the case of asymmetric nuclear 



matter 
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Phase transitions are related to the thermodynamical instabili- 



ties that a system can present. It has been usually argued that asymmetric nuclear matter 
presents two independent types of instabilities: a mechanical (or isoscalar) instability associ- 
ated to density fluctuations which conserves the proton fraction and a chemical (or isovector) 
instability, related to fluctuations in the proton fraction, occurring at constant density. Nev- 



ertheless, it has recently been shown 0,0 



[l3l | that actually asymmetric nuclear matter 
presents only one type of instability whose direction is dominated by total density fluctua- 
tions which lead to a liquid-gas phase separation with restoration of the isospin symmetry 
in the liquid dense phase. This is known as the isospin distillation or fragmentation effect 
[l4 ]. where large droplets of high density symmetric matter are formed in a background of 
a neutron gas with a small fraction of protons. 

There have been systematic analysis of the stability conditions of asymmetric nuclear 
matter against the liquid-gas phase transition by using different approaches, mean field 
calculations with effective forces of the Skyrme or Gogny type 12, 15, Q, l^j], relativistic 



mean field calculations using constant and density dependent coupling parameters [18l.ll9l.l20. 



21 



221 ]. or very recently, a Dirac-Brueckner-Hartree-Fock approach with realistic nucleon- 



nucleon potentials 



]. In the present work, we perform a similar study starting from 
a realistic interaction, namely the Argonne V18 24j plus a three body force of Urbana 
type, in the framework of the Brueckner-Hartree-Fock (BHF) approximation. After a brief 
description of the BHF approach, we will summarize the stability conditions and, then, 



discuss and compare our results with those obtained with the Skyrme force SLy230a 
and the relativistic mean field models NL3 26| and TW [271 ] . 
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Our microscopic many-body scheme starts with the construction of the G-matrix, which 
describes the effective interaction between two nucleons in the presence of a surrounding 
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medium. It is obtained by solving the well known Bethe-Goldstone equation 

(kin; k 2 r 2 | G(u) | k 3 r 3 ; k A n) = (kin; k 2 r 2 \ V \ k 3 T 3 ; k 4 n) 

+ y^ikm; k 2 r 2 | V | hn; kfo) ® TlT J — —(hn; kjTj \ G(u) I k 3 r 3 ; k 4 n) , (1) 

*-r? oj - E n - E T . + irj 

where, r = n,p, indicates the isospin projection of the two nucleons in the initial, interme- 
diate and final states, k are their respective linear momenta, V denotes the bare interaction, 
Q TiT j is the Pauli operator which allows only intermediate states compatible with the Pauli 
principle, and u, the so-called starting energy, corresponds to the sum of non-relativistic 
energies of the interacting nucleons. The single-particle energy E T of a nucleon with mo- 
mentum k is given by 

E T (k) = ^ + Re[U T (k)} , (2) 

where the single-particle potential U T (k) represents the mean field "felt" by a nucleon due 
to its interaction with the other nucleons of the medium. In the BHF approximation, U(k) 
is calculated through the "on-shell energy" G-matrix, and is given by 

U T (k) = ^n Tl (k')(kr; fcV | G{uj = E T {k) + E T t(k')) \ kr; k'r') A (3) 

r'=n,p k' 

where 

( 1, if k < k F 

n T (k) = < 

I 0, otherwise 

is the corresponding occupation number, and the matrix elements are properly antisym- 
metrized. We note here that the so-called continuous prescription has been adopted for 
;he sing le-particle potential when solving the Bethe-Goldstone equation. As shown in Refs. 



28 



291 ]. the contribution to the energy per particle from three-hole line diagrams is dimin- 



ished in this prescription. 
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All the BHF calculations carried out in this Letter use the realistic Argonne V18 
nucleon-nucleon interaction supplemented with a three-body force of Urbana type which 
for the use in BHF calculations was reduced to a two-body density dependent force by 
averaging over the third nucleon in the medium [30[. This three-body force contains two 
parameters that are fixed by requiring that the BHF calculation reproduces the energy and 
saturation density of symmetric nuclear matter, see also Ref. jsij for a recent analysis of the 
use of three-body forces in nuclear and neutron matter. The energy density, E/V, is easily 
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calculated once a self-consistent solution of Eqs. (P - [3]) is achieved: 

r^^E^H' (4) 

T=n,p k ' 

The many-body problem at finite temperature has been considered by several authors 



within different approaches, such as the finite temperature Green's function method 32, 



33 



34j . l35l |. the thermal field method 



Brueckner-Bethe-Goldstone 
limit. Baldo and Ferreira 



or the Bloch-De Dominicis (BD) diagrammatric 
expansion [37]. The latter, represents the "natural" extension to finite temperature of the 

BBG) expansion which is recovered in the zero temperature 
301 ] showed that the dominant terms in the BD expansion were 
those that correspond to the zero temperature of the BBG diagrams, where the temperature 
is introduced only through the Fermi-Dirac distribution 

fr{K T) = l + exp([E T (k,T)-^ T (T)}/T) ' (5) 
Therefore, at the BHF level, finite temperature effects can be introduced in a very good 
approximation just replacing in the Bethe-Goldstone equation: (i) the zero temperature 
Pauli operator Q TT i = (1— n T )(l— n T i) by the corresponding finite temperature one Q TT i(T) = 
(1 — / T )(l — f T i), and (ii) the single-particle energies E T {k) by the ones obtained when n T (k) is 
replaced by f T (k, T) in Eqs. ([2]) and (J3]) • In this case, the self-consistent process implies that 
together with the Bethe-Goldstone equation and the single-particle potential, the chemical 
potential Hr(T) of the nucleon, must be extracted at each step of the iterative process, since 
both the Bethe-Goldstone equation and the single-particle potential depend implicitly on 
it. 

Once a self-consistent solution is obtained the free-energy density is determined by the 
usual thermodynamical relation, 

where E /V is evaluated from Eq. (BJ replacing n T (k) by f T (k,T), and the entropy density 
is calculated in the quasi-particle approximation 



S 



"77 E E^' T )HU{k, T)) + (1 - f T (k, T))ln{\ - f T (k, T))] . (7) 
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We have now all the necessary ingredients to analyze the stability of asymmetric nuclear 
matter against phase separation. The stability of the system is guaranteed if the free- 
energy of a single phase is lower than the free-energy in any two-phases configurations. This 



condition is fulfilled if the free-energy density is a convex 



densities p n and p p , that is, if the curvature matrix 



hnction of the neutron and proton 



hi 



C - = (|^-) T ' •'''' = n ' p (8) 

is positive-definite. This demands that both its trace and its determinant should be positive, 
i.e., 

Tr{Cij) = A + + A_ > 0, Det(Cij) = A+A_ > , (9) 

where 

X± = \ ( Tr(C<i) ± ^ Tr ^) 2 ~ ADet ^)) ( 10 ) 
are the two real eigenvalues of the curvature matrix associated to the eigenvectors (5p„ , Spp) 



5pf _ A± - C 
Spf Cji 
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i,j=p,n. (11) 



Stability requires, therefore, that both eigenvalues should be positive. In asymmetric nuclear 
matter A + is always positive, and only A_ can eventually become negative, signaling the 
beginning of the instability and the phase separation [121 ] . It turns out that the magnitude 
of A + exceeds always that of A_ (A + > |A_|) and therefore the trace appears to be always 
positive. This is due, as it has been already pointed out by Margueron and Chomaz in Ref . 



12j , to the fact that: (i) the dominant dependence of the energy on the isospin asymmetry, 
P = (p n — p p )/ (p n + p p ), is essentially quadratic, and (ii) the symmetry energy is a positive 
increasing function of the total density, p = p n + p p . Consequently, the spinodal region at 
a given temperature will be just determined by the values of neutron and proton densities 
which make the determinant of the curvature matrix negative. 

Before presenting our results, we will give some technical details concerning the evaluation 
of the second derivatives of the free-energy density needed to construct the curvature matrix. 
As it has already been pointed out, the dominant dependence of the energy on the isospin 
asymmetry is essentially quadratic, therefore, at a given temperature, one can separate 
in good approximation the density and isospin asymmetry dependence of the free-energy 
density assuming that 

F(p,P)^F (p)+F sym (p)P 2 , (12) 

where J-'o(p) = JF(p, 0) is the free-energy density of symmetric nuclear matter, and T sym = 
jF(p, 1) — jF(p, 0). Now, expressing all the derivatives with respect to neutron and proton 



densities in terms of derivatives with respect to p and /3, and assuming Eq. (fT2"I) . one can 
perform all the derivatives with respect to (3 analytically and write the matrix elements of 
the curvature matrix as a function only of J-q(p), F symi and their first and second order 
derivatives with respect to p, 

C nn = r'o(p) + r> sym (p)P 2 + mi ~ P) F'sy m (p) + ~ ~ ' ^ ~ ~ ^ 'sym(p) , (13) 

Cnp — Cpn — F"o{p) + sym (p) [3 2 sym{p) H T sym{p) > (14) 

C pp = T o(p) + F sym{p)p -F symKP) H ^ Fsym{P) ■ (15) 

The derivatives J 7 " o(p) , J 7 ' S ym(p) and T" sym {p) can also be done analytically if one assumes 
a particular functional form of the free-energy density in symmetric and neutron matter. In 
our case, we have chosen to fit the quantity p d(F / A) / dp for symmetric and neutron matter 
and to express all the derivatives of J-"o(p) and T sym {p) with respect to p in terms of this 
quantity and its derivatives. The reason for such a choice is double, first p d(F / A) / dp goes 
in both cases, symmetric and neutron matter, to the temperature in the limit of very low 
densities, so one has a fix point in the fit, and second, this choice turns out to be numerically 
very stable. 

In the next we present the results of our microscopic Brueckner approach, and compare 
them with those obtained with the Skyrme force SLy230a and the relativistic mean field 
models NL3 and TW. We start by showing in Fig. [1] symmetric nuclear matter isotherms 
obtained in our BHF approach with Argonne VI 8 plus a Urbana type three body force 
for several temperatures ranging from T = up to the critical temperature T c = 17.5 
MeV (indicated by the solid dot). The border of the spinodal zone is delineated by the 
black dashed curve. As it has been already said, the spinodal region is determined by the 
values of p n and p p for which Det(Cy) < 0. In symmetric nuclear matter this region is 
exactly the one where the pressure decreases with density, and therefore where the system 
is mechanically unstable. As a general feature, we can observe that the spinodal instability 
starts at zero density, temperature and pressure, and disappears well before the saturation 
density p = 0.182 fm -3 is reached. 

Let us consider now the most general case of asymmetric nuclear matter. Fig. [2] depicts 
the projection of the spinodal contours in the neutron-proton density plane (left panel) and 
in the density-/? plane (right panel) for several temperatures ranging from T = up to 

6 



T = 17 MeV. The qualitative behaviour is very similar to the one observed in previous 



ca 
16 



culations using different relativistic or non relativistic mean field models 
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17 
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H H, 



c 



(16) 



22l | . As the temperature increases, the instability region, defined 
by the inner part of the spinodal boundaries, shrinks up to the point corresponding to 
the critical temperature of symmetric nuclear matter T c . Note that the isospin symmetry 
of the nucleon-nucleon interaction ensures that the contours are symmetric respect to the 
interchange of p n and p p (left panel) or the interchange of (3 and —(3 (right panel). Notice 
also (see right panel) that for a given temperature, the total density at which the spinodal 
instability appears decreases when the isospin asymmetry increases, similar to what happens 
with some Skyrme (SLy230a, SGII) and Gogny (DIP) forces (12J. We have also considered 
the presence of a degenerated ultrarelativistic free Fermi gas of electrons in the system with 
a density equal to that of the protons in order to keep charge neutrality. In this case, the 
curvature matrix C becomes 

' d 2 Tldp 2 n d 2 F/d Pn d Pp 
d 2 F/dp p dp n d 2 T/dp 2 p + d/i e /dp e 
where p e is the electron chemical potential. We have found that in the density region of 
interest for our study, the electron contribution to the curvature matrix, dp e /dp e , is so large 
that it leads to a complete suppression of the instability region already at zero temperature. 
Nevertheless, this result is model dependent. For example, in Ref. 17| . using Skyrme forces 
SLy230a, SGII and SIII, it was also found a complete suppression of the instability region, 
while for instance in Ref. 20], where the relativistic mean field model NL3 was used, it 
is reported a strong reduction of the instability region in the case of nuclear matter with 
electrons, but not its complete suppression. 

Let us compare now our Brueckner prediction for the spinodal region with those obtained 
with the Skyrme force SLy230a and the relativistic mean field models NL3 and TW. This 
comparison is presented in Fig. [3] for a fixed temperature of 10 MeV. The results for the 
Skyrme force SLy230a have been calculated also by us, whereas those of the relativistic 



2l| . Clearly the larger spinodal 



mean field models NL3 and TW have been taken from Ref. 
region, which almost envolves all the others, corresponds to the Brueckner one. This can be 
easily understood noting that in the symmetric nuclear matter case the larger the critical 
temperature and saturation density at zero temperature are, the larger the spinodal region is. 
Although in asymmetric nuclear matter the onset of the instability region will be sensitive to 



the density dependence of the symmetry energy and its derivatives, one can still expect that 
the model that predicts the larger critical temperature and saturation density for symmetric 
matter will predict the larger spinodal instability region in asymmetric matter. In fact, our 
Brueckner calculation predicts the larger critical temperature, T c = 17.5 MeV and saturation 
density, po — 0.182 fm~ 3 . The closest curve to the Brueckner prediction corresponds to 
the one obtained for the non-relativistic Skyrme force SLy230a, which predicts a critical 
temperature of T c = 14.6 MeV and a saturation density po = 0.160 fm -3 . Finally, the critical 
temperatures for NL3 and TW are 14.5 and 15.1 MeV respectively, while the saturation 
densities are quite similar, 0.148 fm -3 for NL3 and 0.153 fm -3 for TW, predicting smaller 
spinodal regions. 

Usually, it has been argued that asymmetric nuclear matter presents two independent 
types of instabilities: a mechanical (or isoscalar) instability conserving the proton fraction 
and a chemical (or isovector) instability occurring at constant density. Nevertheless, recently 
Margueron and Chomaz 12] (see also Refs. [lo| and 13]) have shown that in fact instabilities 
in asymmetric nuclear matter appear as a mixture of density and proton fraction fluctuations, 
and, therefore, they cannot be separated into mechanical or chemical instabilities: chemical 
(mechanical) instabilities do not only involve changes in the proton fraction (density) but 
also changes in the density (proton fraction). There is only one type of instability whose 
direction is given by the eigenvector (Sp~,Sp~) associated with the negative eigenvalue A_ 



[lo, 
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131 ] . If 5p p /5p n = p p / p n then the instability preserves the ratio between protons and 



neutrons at which the system was prepared, and its nature is purely mechanical, whereas 
if Spp = —Sp~ then the total density remains constant and therefore the instability is 
purely chemical. In general, however, the nature of the instability will never be neither 
purely mechanical nor purely chemical, but it will appear as a mixture of both which will 
)e predominantly of isoscalar type if 5p~/5p~ > or of isovector type if 5p~/5p~ < 



m Ha, iii 

Figure |H shows the ratio Sp~/5p~ as a function of the proton fraction (upper panels) 
for several temperatures ranging from T = 0uptoT=17 MeV and three fixed densities: 
p = 0.03 (left panel), 0.06 (middle panel) and 0.09 (right panel) fm -3 . The same ratio is 
showed in the lower panels as a function of the total density for the same set of temperatures 
and proton fractions: x p = 0.05 (left panel), x p = 0.2 (middle panel) and x p = 0.3 (right 
panel). Note that in all cases Sp~/Sp~ is positive, indicating that the instability is always 



dominated by the total density fluctuations (i.e., it is of isoscalar type) even for large 
asymmetries, in ag — with thefts of MS . Q, H and Q. Th 18 is expected 
because the symmetry energy in our approach is always a positive increasing function of 
density, and therefore we do not expect the instabilities to be dominated by proton fraction 
fluctuations. In particular, we notice that for symmetric matter (x p = 0.5) 5p~/5p~ = 1, 
meaning that in this case the instability occurs in the pure isoscalar direction and that 
matter behaves as a one component system. Moreover, it can be observed that 5p~ /5p~ is 
larger than the ratio between protons and neutrons p p / p n (see upper panels), which tells 
us that the instability drives the dense phase (liquid) of the system towards a more isospin 
symmetric region of the p n -p p plane. The conservation of the total particle number enforces 
the light phase (gas) to become more neutron rich leading to the so-called isospin distillation 
or fragmentation effect [fjj]. The ratio 5p p /5p~ can be then considered a measurement 
of the efficiency of isospin symmetry restoration in the liquid phase: the larger its value 
the greater the efficiency In general, we observe that whereas the restoration of isospin 
symmetry becomes more and more efficient when the proton fraction increases (5p~/5p~ 
increases), it becomes less efficient when the temperature and the density increase (5p~ /6p~ 
decreases). Only for very low densities (p < 0.02 fm -3 ) this efficiency grows with increasing 
density. 

Finally, we compare in Fig. [5] our Brueckner prediction for the isospin symmetry restora- 
tion with those of the Skyrme SLy230a force and the relativistic mean field models NL3 and 
TW. In general, one can say that in the Brueckner approach the restoration of isospin sym- 
metry is less efficient (5p~/5p~ is smaller) than in the case of Skyrme or relativistic mean 
field, except for very low proton fractions. We also note (see lower panels) that whereas 
the Brueckner approach, the Skyrme SLy230a and the TW model share the same qualita- 
tive behaviour (Sp~/Sp~ decreases with density), the NL3 model predicts the opposite one 
(6p~/5p~ increases with density). As it has been discussed by Avancini et al. in Ref. 211 ] . 
this is essentially due to the presence of the rearrangement contribution in the evaluation 
of the neutron and proton chemical potentials. This contribution is present in the BHF 
approach, in the Skyrme force and in the TW model, but it is completely absent in the NL3. 
In fact, the authors of Ref. 21[ showed that by removing this contribution from the TW 



model, the results predicted by this model (TW with no rearrangement) show a behaviour 
similar to NL3. 
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To enlight the analysis of the distillation effect, we report in Fig.[6]the density dependence 
of the symmetry energy predicted by the different models. Nevertheless, it is clear from 
Eqs. (I13II14|I15I) that the matrix elements of the curvature matrix, and therefore the ratio 
Sp~/5p~, depend not only on the value of the symmetry energy but also on its first and 
second derivatives, and therefore to draw a conclusion on the distillation effect for the 
different effective interactions is not straightforward and would require a more sophisticated 
and careful analysis, than just considering the value of the symmetry energy. This analysis 
is beyond the scope of the present work and it will be addressed in the future. The figure 
also serves to illustrate that whereas the BHF approach, the Skyrme SLy230a force and the 
TW model predicts a qualitatively similar behaviour of the symmetry energy as a function 
of density, the NL3 model clearly shows a different one, being almost linear. This different 
behaviour of the symmetry energy is in agreement with the fact, already discussed, that 
NL3 model predicts also a behaviour of the ratio Sp~/5p~ qualitatively different from the 
one predicted by the other three models. 

Summarizing, we have studied the spinodal instabilities of asymmetric nuclear matter 
at finite temperature within the microscopic Brueckner-Hartree-Fock approach using the 
realistic Argonne V18 potential plus a three-body force of Urbana type which guarantees 
the reproduction of the saturation density in symmetric nuclear matter. The results are 
qualitatively similar to previous analysis performed with non-relativistic and relativistic 
effective interactions. In particular, we have compared our results with those obtained with 
the Skyrme force SLy230a and the relativistic mean field models NL3 and TW, finding 
that the Brueckner approach predicts a larger spinodal region. This is mainly related to 
the fact that our Brueckner calculation with Argonne V18 plus a Urbana type three-body 
force predicts values for the critical temperature and the saturation density of symmetric 
nuclear matter larger than the ones provided by those effective interactions. We have found 
that the instability is always dominated by total density fluctuations, in agreement with 
previous results of other authors 1(], 12, 21]. Finally, we have also studied the restoration 



of the isospin symmetry in the liquid phase, finding that whereas it becomes more and more 
efficient when proton fraction increases, it becomes less and less efficient with the increase of 
the temperature and the density. In general, we found that the BHF results are comparable 
to those obtained with the Skyrme and the relativistic mean field models, although the 
restoration of isospin symmetry is not so efficient in this case. 
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FIG. 1: (Color online) Symmetric nuclear matter isotherms for several temperatures. The border 
of the spinodal zone is delineated by the black dashed curve; uniform matter below this curve is 
mechanically unstable. The critical point, indicated by the solid dot, corresponds to a temperature 
of T c = 17.5 MeV. 
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FIG. 2: (Color online) Projection of the spinodal contour in the neutron-proton density plane (left 
panel) and in the density-/? plane (right panel) for temperatures ranging from T = up to T = 17 
MeV. 
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FIG. 3: (Color online) Comparison of the spinodal contour at T = 10 MeV predicted by our BHF 
approach (black solid line) with those obtaine d by us with the Skyrme force SLy230a (red dotted 
line) and those obtained by Avancini et al. [2l| using the relativistic meand field models NL3 
(green dashed line) and TW (blue dot-dashed line). 
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FIG. 4: (Color online) Upper panels: Ratio 5p~ /5p~ as a function of the proton fraction for several 
temperatures and fixed densities, p = 0.03 (left panel), 0.06 (middle panel) and 0.09 (right panel) 
fm~ 3 . Black dotted line shows the ratio between protons and neutrons p p /p n - Lower panels: Same 
ratio as a function of the total density for the same set of temperatures and fixed proton fractions, 
x p = 0.05 (left panel), x p = 0.2 (middle panel) and x p = 0.3 (right panel). 
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FIG. 5: (Color online) Comparison of the ratio 5p p /5p n at T = predicted by our BHF approach 
(black solid lines) with those obtained by us with the Skyrme force SLy230a (red dotted lines) and 



those obtained by Avancini et al. 



2l| using the relativistic mean field models NL3 (green dashed 



lines) and TW (blue dot-dashed lines). As in the previous figure, in the upper panels this ratio is 
plotted as a function of the proton fraction for three fixed densities, p = 0.03,0.06 and 0.09 fm~ 3 , 
whereas its density dependence is shown in the lower panels for fixed proton fractions, x p = 0.05 
(left panel), x p = 0.2 (middle panel) and x p = 0.3 (right panel). 
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FIG. 6: (Color online) Comparison of the symmetry energy as a function of the density obtained 
within our BHF approach (black solid line), with those obtained with the Skyrme force SLy230 
(red dotted line) and the relativistic mean field models NL3 (green dashed line) and TW (blue 
dot-dashed line). 
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